01 引言

上一篇推文【Python量化基础】时间序列的自相关性与平稳性着重介绍了时间序列的一些基础概念，包括自相关性、偏自相关性、白噪声和平稳性，以及Python的简单实现。本文在此基础上，以沪深300指数收益率数据为例，探讨如何使用Python对平稳时间序列进行建模和预测分析。时间序列经典模型主要有自回归模型AR，移动回归模型MA，移动自回归模型ARMA，以及差分移动自回归模型ARIMA，今天主要介绍这四种模型的基本原理以及Python的实现步骤。

为了阅读体验效果更好，请阅读公众号原文：https://mp.weixin.qq.com/s/vLuvBXpF96HzPX6mdf7x8w

02 AR模型

AR模型全称为Autoregressive Models，即自回归模型，用于刻画因变量能由它的多个滞后项表示。p阶自回归模型可以写成：

下面模拟一个AR(1)模型。

importpandasaspd

importnumpyasnp

importstatsmodels.tsa.apiassmt

#tsa为Time Series analysis缩写

importstatsmodels.apiassm

importscipy.statsasscs

fromarchimportarch_model

#画图

importmatplotlib.pyplotasplt

importmatplotlibasmpl

%matplotlib inline

#正常显示画图时出现的中文和负号

frompylabimportmpl

mpl.rcParams['font.sans-serif']=['SimHei']

mpl.rcParams['axes.unicode_minus']=False

#先定义一个画图函数，后面都会用到

defts_plot(data, lags=None,title=''):

ifnotisinstance(data, pd.Series):

data = pd.Series(data)

#matplotlib官方提供了五种不同的图形风格，

#包括bmh、ggplot、dark_background、fivethirtyeight和grayscale

withplt.style.context('ggplot'):

fig = plt.figure(figsize=(10,8))

layout = (3,2)

ts_ax = plt.subplot2grid(layout, (0,0), colspan=2)

acf_ax = plt.subplot2grid(layout, (1,0))

pacf_ax = plt.subplot2grid(layout, (1,1))

qq_ax = plt.subplot2grid(layout, (2,0))

pp_ax = plt.subplot2grid(layout, (2,1))

data.plot(ax=ts_ax)

ts_ax.set_title(title+'时序图')

smt.graphics.plot_acf(data, lags=lags, ax=acf_ax, alpha=0.5)

acf_ax.set_title('自相关系数')

smt.graphics.plot_pacf(data, lags=lags, ax=pacf_ax, alpha=0.5)

pacf_ax.set_title('偏自相关系数')

sm.qqplot(data, line='s', ax=qq_ax)

qq_ax.set_title('QQ 图')

scs.probplot(data, sparams=(data.mean(),

data.std()), plot=pp_ax)

pp_ax.set_title('PP 图')

plt.tight_layout()

return

# 模拟AR(1) 过程 

#设置随机种子（括号里数字无意义）

np.random.seed(1)

#模拟次数

n=5000

#AR模型的参数

a =0.8

#扰动项为正态分布

x = w = np.random.normal(size=n)

fortinrange(1,n):

x[t] = a*x[t-1] + w[t]

#画图

ts_plot(x, lags=30)

模拟的AR(1)模型是正态的。自相关系数图（ACF）显示滞后值之间存在显著的序列相关性，偏自相关系数图（PACF）则显示在滞后1期时截尾（迅速降为0）。下面使用statsmodels构建AR(p)模型，先用AR模型拟合上述模拟的数据，并返回估计的系数参数），然后选择最佳滞后阶数，最后与原模型设置对比看是否选择了正确的滞后项。假如AR模型是正确的，那估计的系数参数将很接近真实的系数0.8，选择的阶数也会等于1。

#估计数据的AR模型参数和滞后阶数

defsimu_ar(data,a,maxlag=30,true_order =1):

'''data:要拟合的数据；a为参数,可以为列表；maxlag:最大滞后阶数'''

# 拟合AR(p)模型

result = smt.AR(data).fit(maxlag=maxlag, ic='aic', trend='nc')

#选择滞后阶数

est_order = smt.AR(data).select_order(maxlag=maxlag,

ic='aic', trend='nc')

#参数选择标准ic : 有四个选择 {‘aic’,’bic’,’hqic’,’t-stat’}

#趋势项：trend：c是指包含常数项，nc为不含常数项

#打印结果

print(f'参数估计值：{result.params.round(2)}，

估计的滞后阶数：{est_order}'

)

print(f'真实参数值：{a}，真实滞后阶数{true_order}')

simu_ar(x,a=0.8)

参数估计值：[0.8]，估计的滞后阶数：1

真实参数值：0.8，真实滞后阶数 1

看下如何用AR(p)模型来拟合沪深300的对数收益

# Select best lag order for hs300 returns

importtushareasts

token='输入token'

pro=ts.pro_api(token)

df=pro.index_daily(ts_code='000300.SH')

df.index=pd.to_datetime(df.trade_date)

deldf.index.name

df=df.sort_index()

df['ret']=np.log(df.close/df.close.shift(1))

max_lag =30

Y=df.ret.dropna()

ts_plot(Y,lags=max_lag,title='沪深300')

result = smt.AR(Y.values).fit(maxlag=max_lag, ic='aic', trend='nc')

est_order = smt.AR(Y.values).select_order(maxlag=max_lag,

ic='aic', trend='nc')

print(f'沪深300拟合AR模型的参数：{result.params.round(2)}')

print(f'沪深300拟合AR模型的最佳滞后阶数{est_order}')

沪深300拟合AR模型的参数：[ 0.03 -0.03  ...]

沪深300拟合AR模型的最佳滞后阶数 15

最好的阶数选择是15或者说有15个参数！任何模型有这么多参数在实际中不可能有用。显然有比这个模型更好的模型可以解释沪深300收益率走势。

03 MA模型

MA(q)模型与AR(p)模型非常相似。不同之处在于，MA(q)模型是对过去的白噪声误差项的线性组合，而不是过去观测的线性组合。MA模型的动机是我们可以直接通过拟合误差项的模型来观察误差过程中的“冲击”。在一个AR(p)模型中，通过在一系列过去的观察中使用ACF间接观察到这些冲击。MA(q)模型的公式是：

下面使用Python模拟MA(1) 过程。

#这里使用arma模型进行模拟，设定ar阶数为0，即得到ma模型

alphas = np.array([0.])

betas = np.array([0.6])

ar = np.r_[1, -alphas]

ma = np.r_[1, betas]

#模拟MA的样本数据

ma_sample = smt.arma_generate_sample(ar=ar, ma=ma, nsample=1000)

ts_plot(ma_sample, lags=30,title='MA(1)模型')

ACF函数显示滞后1阶系数显著异于0，表明MA（1）模型适合拟合的数据。

# 对上述模拟数据进行ARMA模型拟合

max_lag =30

result = smt.ARMA(ma1, order=(0,1)).fit(maxlag=max_lag,

method='mle', trend='nc')

print(result.summary())

模型估计d 滞后系数为0.6277，与真实值0.6比较接近。注意到，95％置信区间确实包含该真实值。

下面尝试用MA(3)模型去拟合沪深300股价的对数收益，但这次并不知道真实的参数值。结果显示，拟合的残差自相关系数和偏自相关系数比较符合白噪声过程，但由于存在厚尾，MA模型并不是预测沪深300未来回报的最佳模型。

max_lag =30

result=smt.ARMA(Y.values,order(0,3)).fit(maxlag=max_lag,

method='mle', trend='nc')

print(result.summary())

resid=pd.Series(result.resid,index=Y.index)

ts_plot(resid, lags=max_lag,title='沪深300指数MA拟合残差')

04 ARMA模型

ARMA模型全称为自回归移动平均模型Autoregressive Moving Average Models - ARMA(p, q)，是AR（p）和MA（q）模型之间的结合，从金融的角度理解，AR和MA模型的理论意义在于：AR（p）模型试图捕捉（解释）交易市场中经常观察到的动量和均值回复效应。MA（q）模型尝试捕捉（解释）在白噪声条件下观察到的冲击效应。这些冲击效应可以被认为是影响观察过程的意外事件。ARMA模型的弱点在于忽视了大多数金融时间序列中的波动聚集效应。模型的公式可以表示为:

print(result.summary())

# 下面使用ARMA(2, 2) 模型进行模拟分析

max_lag =30

n =5000

burn = int(n/10)

alphas = np.array([0.5,-0.25])

betas = np.array([0.5,-0.3])

#注意ar模型1代表0阶（自身），然后在其他系数前加负号

ar = np.r_[1, -alphas]

ma = np.r_[1, betas]

arma22 = smt.arma_generate_sample(ar=ar, ma=ma, nsample=n, burnin=burn)

_ = ts_plot(arma22, lags=max_lag)

result = smt.ARMA(arma22, order=(2,2)).fit(maxlag=max_lag,

method='mle', trend='nc', burnin=burn)

结果显示模型估计的参数与真实参数基本上吻合。下面使用ARMA模型来拟合沪深300的收益数据。ACF和PACF没有显示出明显的自相关性。QQ和概率图显示残差大致为正态分布但厚尾。总体而言，这个模型的残差看起来不像白噪声，说明模型还是没有很好的拟合其波动性特性。

#不事先确定滞后阶数，而是通过信息准则选择最佳的滞后阶数

#先将初始值设置为无穷大

best_aic = np.inf

best_order =None

best_mdl =None

rng = range(5)

foriinrng:

forjinrng:

try:

tmp_mdl = smt.ARMA(Y.values, order=(i,j))

.fit(method='mle', trend='nc')

tmp_aic = tmp_mdl.aic

iftmp_aic